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INELASTIC BUCKLING OF NON- UNIFORM COLUMNS 


John E. Goldberg, ! M. ASCE, John L. Bogdanoff ,2 and Hsu Lo® 
(Proc. Paper 943) 


SYNOPSIS 


A numerical procedure is presented for calculating the critical loads of 
non-uniform pin-ended columns in the elastic or inelastic range. The proce- 
dure is based upon transformation of the differential equation to an integral 
equation which is solved numerically for trial values of the critical load. 


If an initially straight pin-ended column is subjected to a slowly increasing 
axial load, the column will buckle when the load exceeds the value of P cor- 
responding to the first eigenvalue of the following differential equation togeth- 
er with appropriate boundary conditions: 


4 
+ py = 0 


(1) 
2 
ax’ t 
where E, = E, (x5P) = tangent modulus 


I = I(x) = moment of inertia of cross section 
P = axial load 
x = longitucinal coordinate 
y = leteral deflection 
Also L = length of colum,. 
Note: Discussion open until September 1, 1956. Paper 943 is part of the copyrighted 
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The above equation merely states that in some particular configuration, 
y = y(x), which is arbitrarily small, the column is in a state of equilibrium 
when the applied load has been increased slowly to the critical load. Taking 
the internal bending moment as 


2 

(2) 
. 
ax 


and equating this expression to the applied moment, P times y, leads to (1). 
For the case of the pin-ended column, the boundary conditions may be 
taken as 


y(0) = 0 (3a) 


y(L) = 0 (3b) 
or, alternatively, since the bending moment vanishes at the ends of the column, 
(3c) 


dx” 


9 (3d) 


It is convenient to rewrite the above equations in terms of a non- 
dimensional independent variable 


and it is easily seen that 


Consequently, the use of (4) and (4a) in (1) results in the equation, 
2 


“3 77 


with boundary conditions 


4 
ee 
4 s = (4) 
2 
5 
2 
ds t 
y(0) = 0 (5a) 


ASCE GOLDBERG, et al. 


or, alternatively, 


(5c) 
(5d) 


Equation (5) is readily transformed to an integral equation by the following 
device. We substitute for the second derivative of y a new dependent variable 
Y(s), so that 


(6a) 


ds 


y =f, - (6c) 


We observe that the lact of (6) and the boundary condition (5a) at s = 0 require 
that 
Cc, = 0 


2 


The substitution of (6) into (5) therefore yields the following Volterra integral 
equation of the second kind, 


2 
¥(s) = Et + + (7) 


while the boundary condition at s = 1 becomes 


¥(1) = 0 (Ta) 


By (7) it is found that Y(0) = 0 which, of course, agrees with (5c). 

A solution of (7), linear in Cj, can be found by numerical methods de- 
scribed below for any value cf P. The values of P associated with solutions 
satisfying the boundary condition, (7a), represent estimates of critical loads 
corresponding to the various modes of buckling of the column. The lowest of 
these P’s greater than zero is the critical load of the column. 

In terms of C, and for a given value of P, (7) represents an initial value 
problem for which an approximate solution may be found by the following 
method. We divide the range of integration (0,1) into n equal segments with 
boundary and division points located at sp, 53, 89,......, S,_1>Sp- The 
points sg and s, correspond to the ends of the column. The le of each 
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segment is b = 1/n. At the j-th division point Y(s) takes on the approximate 
value Y;. Between division points, an interpolation of Y(&) is used, thus per- 
mitting evaluation of the sectionalized integrals of (7). If the column is divid- 
ed into a sufficient number of segments, the assumption that Y(&) varies 
linearly over each segment will give good results. 

Assuming linear variation over each segment, the approximate form of 
Y(€) for integration in (7) takes on the polygonal character indicated in Fig. 1. 
With this assumption, it is easily seen that 


Jj 
ft Y(E) dk, = + bY, 


We » Yor 


Ky + (d- @) 


Substituting (8) into (7), collecting terms in Yj, and making use of the fact 
that Yo = 0 leads to 


1 


r 2,2 2 2 j-1 j-1 
- Fu 
which may be written as 


If the summations in (9) are transposed to the left-hand side, it is seen that 
the solution must be linear in C;. Consequently, for the present purposes, it 
suffices to select any convenient value for C;, say unity. 


Suggested Procedure 


The following procedure is suggested for the determination of the critical 
load of a non-uniform column. The procedure is essentially a trial-and-error 
method. 


1) Divide the column into n equal segments and tabulate area, A, and mo- 
ment of inertia, I, at each division point and at the terminal point. 


2) Select a trial value of the axial load. 


3) Compute the average stress at each division point and at the terminal 
point, and determine the tangent modulus for each point at the applied 
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stress. The tangent modulus may be determined from a tangent modu- 
lus curve or from a stress-strain curve, or it may be determined by 
means of an analytical expression for the stress-strain curve. 


4) Using the following formula, 


+ Sux 
Bix 
j 2. 


n E,I 


Y (9a) 


compute Yj at j=1,2,...,m. This may be done by desk computing 
methods or by digital computer. 


5) If Y,, is found to be zero, then the trial load which has been used is a 
critical load. It is the critical load corresponding to the lowest mode 
of buckling if Yj has the same algebraic sign throughout. 


6) If Y, is not zero (or approximately so), a new trial value of P should 
be selected and steps 2 to 5 repeated. If Y changes sign before the end 


point, n, is reached, this implies that a point of inflection or a hinge Hi 
exists within the column and, therefore, the previous trial load is the 

critical load for a shorter column. Consequently, for the next trial, a 

lesser value of P should be used. Conversely, if Y does not change sign 

within the column, this corresponds to buckling of a somewhat longer " 


column and, consequently, a larger value of P should be selected for the 
next trial. 


It is, of course, unlikely that the boundary condition at s, will be satisfied 
by the first trial. Since a second trial will be necessary, it is desirable to use 
the results of the previous calculation as a guide in selecting a new trial value 
of the axial load. By interpolation or extrapolation, the approximate point at 
which Y becomes zero may be found. If the point at which Y becomes zero is 
close to the end of the column, an estimate of the correction to the trial value 
of the critical load may be obtained by the following device. Let As be the 
distance of this inflection point from sp, being positive if the inflection point 
is found to lie beyond the end of the column, and negative if it is found to lie 
within the column, Then the correction, AP, to the previous trial value, Pj, 
may be taken as 


AP = 2P,As (10) 


and the value to be used as a second trial may be taken as 


Equation (10) is obtained by differentiating the Euler formula with respect 
to the length, assuming that the average tangent modulus will not change 
greatly. Thus 


an = 
4,1 


(12) 
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where the bars over E; and I indicate properties having the character of 
averages, then 


2n°E,T AL 
5 
Li 
or AP = 2P, 4s (13) 


If the additional terms which enter into the total differential of (12) can be 
computed readily, i.e., the terms involving the differentials of the tangent 
modulus and the moment of inertia, a more refined expression for the cor-- 
rection to the previous trial value may be obtained. However, it is doubtful 
whether, under ordinary circumstances, too much effort at this point is jus- 
tified since the final result probably will be obtained by interpolation. Since 
(10) is not in general the total differential of (12), this equation should be used 
only as a means of estimating a correction. The final result is best obtained 
by interpolation or extrapolation from trial results which are reasonably well 
converged. 

A variation of (9) is sometimes convenient for computation. For this pur- 
pose, let the product EI at some station, say sg, be selected as a reference 
value and designated (EjI)g. Also let 


= = r(s,p) 


and let 


The latter represents that ratio of Euler load of a uniform bar having the 
section properties of the reference station to the trial load. Then (9) may be 


written 
2 j-1 
+ = C,jn= j + ky (9b) 
pT 1 1 


Examples and Results 


Example 1. Uniform Bar in Elastic Range 


This example is included merely to establish the potentialities of the 
method with respect to accuracy. The problem was worked using 5, 10, and 
20 intervals. For a uniform bar, rj in (9b) is always equal to unity. As dis- 
cussed above, Cy is taken as unity. 

In each case, a first trial value was taken as five percent greater than the 
Euler load; hence py = 1.05. Second trials were made with po = 0.95. In each 
case, correction formula (13) was also applied to obtain a value which could 
be taken, if desired, for a subsequent trial. Linear interpolation or extrapo- 
lation was used to determine As, realizing, of course, that this is increasingly 
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inaccurate as n decreases and that a second order difference expansion should 
be used for better results. Starting with po for ng and using extrapolation 
based upon second differences to compute As gives a corrected p of 1.018 in- 
stead of 1.0472 as shown for linear extrapolation. 

An alternative procedure which consists merely of interpolating between 
P 1 and pg to find the value of p at which Yp vanishes was also used. These 
interpolated values were taken as a third trial value and the calculations were 
repeated, leading to the results shown in the last table. 


Example 2. Column with Linearly Varying Thickness 


A 1"'-wide solid column having a thickness varying linearly from 1" to 
0.5" and a length of 10" is considered. Values of the tangent modulus are 
taken from Fig. 2. Calculations are made following the scheme given in 
Table II, with n = 10. 

A first trial was made using P = 50,000 lbs. and gave Y,, = -1.603 as com- 
pared with a maximum Y equal to 7.48. This indicates that the trial value of 
P was too large. Using a rough graph of the values of Yg, Yg, and Y, in lieu 
of an interpolation formula, the point of inflection was found to be about one- 
tenth of one segment above the bottom, hence As may be taken as -0.01, and 
by (13), AP = 2(50,000) (-0.01) = -1000 lbs. It is realized, however, that the 
apparent precision is misleading. Nevertheless, the indicated correction was 
applied and the process repeated with P = 49,000 lbs. The calculations are 
shown in Table II. 

It is seen that, for P = 49,000 lbs., Y, = 0.488. Interpolating between this 
result ana the result for P = 50,000 lbs. gives a final critical load of 49,233 
lbs. 


CONCLUSION 


Numerous methods exist for transforming equation (7) into a set of linear 
algebraic equations. In the foregoing, linear interpolation has been used over 
each segment. This has led to a simple formula, (9), for constructing a set 
of algebraic equations whose coefficients form a triangular matrix, and hence 
these equations are easily solved. It may be observed that if a greater num- 
ber of stations is taken in order to improve accuracy, which may be justified 
particularly for irregular columns, the labor increases only proportionately 
and not as a higher power of the order of the matrix. 

In lieu of the linear interpolation which has been employed, it is clearly 
possible to use parabolic, cubic, or higher order interpolation formulas. This 
would require the simultaneous solution, by pairs in the case of the parabolic 
formula, by threes in the case of the cubic formula, and so on, of the final 
equations. An alternative procedure which retains the triangular nature of 
the system of equations and yet provides generally improved accuracy may 
be obtained by using single-step forward-integration based upon any desired 
degree of polynomial. Ultimately, one may construct an interpolating poly- 
nomial of the Lagrange type which covers all division points on the column. 
Such methods tend to lose the computational simplicity of equation (9), and it 
is perhaps preferable to use the linear method with a larger number of sta- 
tions if greater accuracy is desired. 

The method which has been presented may be extended to the case of 
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columns with one or more abrupt changes in section. Different arbitrary 
constants of integration will exist for each portion, but relations between 
them are obtained from the continuity conditions at each step. 
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TABLE I 


Example 1 


FIRST TRIALS, VARIOUS n's, p = 1.05 


Interp.p Y, for 
for Y,2° Interp.p 


-.01193 00016 


201031 «01403 00016 


-.00401 201946 00008 


a 
— = 
216 -.91193 201395 -20231 -.0485 1.0015 
216 -.01031 204104 -.0201 -.0422 1.0078 

216 00401 -0$368 -.0982 -.0172 1.0328 
SECOND TRIALS, VARIOUS n'‘'s, p = 0.95 

10 215 201403 05935 20309 20649 1.0144 
5 215 291946 «10358 20463 20972 1.0472 
ALTERNATE INTERPOLATION PROCEDURE AND KESULTS OF THIRD TRIALS 
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1. POLYGONAL APPROXIMATION FOR Y(s). 


4 
40 60 80 100 120 


ox 107? 
FIG. 2. TANGENT MODULUS CURVE FOR EXAMPLE 2. 
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THE DYNAMIC RESPONSE OF TALL STRUCTURES TO LATERAL LOADS 


L. Schenker,! A.M. ASCE 
(Proc. Paper 944) 


ABSTRACT 


The paper deals with some aspects of the analytical procedures for com- 
puting the response of structures to dynamic forces. It is restricted to the 
consideration of tall structures, such as tier buildings and chimneys, sub- 
jected to lateral forces or displacements. The concept of “stiffness matrix” 
is used to develop a general method for dealing with structures in which ro- 
tation of the joints can take place. The method can be extended to take ac- 
count of the semi-rigidity of connections, the flexibility of foundations and 
deformation of parts of the structure into the inelastic range. The effect of 
damping on the method of analysis is discussed in some detail. 


INTRODUCTION 


A study which has the ultimate objective of increasing the resistance of 
structures to dynamic forces can be divided into four phases. Phase one is 
the determination or the estimate of the systems of forces which it is prac- 
ticable to have structures withstand. This phase must be largely experimen- 
tal and involves the analysis of field data in the light of economic considera- 
tions. 

Phase two is the determination of properties of structures that govern 
dynamic behavior: the distribution of mass, relative stiffness values and the 
nature of damping. 

Phase three concerns the analysis transforming the information from 
phases one and two into structural response in terms of stresses, deflections, 
moments and shears. 

Since the analytical procedures of phase three are bound to be complex and 
time consuming there is a need for simplification leading to procedures suit- 
able for design purposes. This can be considered the fourth phase. 


Note: Discussion open until September 1, 1956. Paper 944 is part of the copyrighted 
Journal of the Engineering Mechanics Division of the American Society of Civil 
Engineers, Vol. 82, No. EM 2, April, 1956. 


1. Member of Technical Staff, Bell Telephone Labs., Inc., Murray Hill, N. J. 
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This paper deals with some aspects of phase three, the analytical proce- 
dures for computing the response of structures to dynamic forces. It is re- 
stricted to the consideration of tall structures, i.e. tier buildings, chimneys 
and towers, subject to lateral forces or displacements. The effect of a num- 
ber of factors, usually neglected in the dynamic analysis of structures, is 
discussed. 


Elastic Structures Without Damping 


The following assumptions are made in this study: 


i) The mass of a structure can be represented by a number of discrete 
masses concentrated at various stations. 


ii) The external forces are considered to be concentrated at the stations. 


It is frequently assumed in the dynamic analysis of tall structures that the 
motion of any one story is resisted only by the spring forces brought into play 
by the relative displacement from the two neighboring stories. This implies 
that the sets of columns between floors act like simple springs connecting the 
masses of the floors, as shown in Fig. 1(a). This assumption is justified only 
if the girders are much more rigid than the columns, so that almost no rota- 
tion takes place at the joints. In the more general case, discussed in this 
paper, when the stiffness of the girders is of the same order as that of the 
columns, the motion of each floor is affected by the motion of every other 
floor. The model in Fig. 1(b) represents this state of affairs. Under this 
condition the equations of motion can no longer be set up by inspection. 

For a given system of external lateral forces F,(t)—(r refers to the sta- 
tions)—the dynamic response of a structure can be calculated, provided that 
the numerical values of all masses My, and the lateral stiffness matrix 


Si) 543 


are known. The lateral stiffness matrix for a structure with n masses is de- 
fined here as a square array of n2 numbers, the r-th row of which represents 
the static forces that have to be applied at the various levels in order to pro- 
duce a unit displacement at the r-th level, but a zero displacement at all other 
levels. 

The importance of the concept of the stiffness matrix [S] is that it is a 
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static property of the structure. Yet, apart from the data on the masses, it 
is all the information that is required to compute the undamped response of 
a structure to any given system of lateral time-varying forces. The refine- 
ment with which [S] is computed or computable governs the accuracy of the 
final answer. For example, if[S] is computed, taking account of the rocking 
at the foundation, lengthening and shortening of columns, rotation of the joints, 
the answer is correspondingly characteristic. Furthermore, it is feasible to 
determine [S] experimentally. 

The equations describing the undamped motion which results from the ap- 
plication of a system of forces F,. are: 


+ 47} (2.2) 


where [M] is the diagonal matrix representing the mass at every station. 
Salvadori(1)* has shown that these simultaneous differential equations can 
be solved conveniently by the transformation 


j=n 
x(t) = .q,(t) (2.3) 
j=l rj j 
where @,.; is the normalized shape of the j-th mode, corresponding to the 
natural reular) frequency w;. The generalized coordinates qj are deter- 


mined from the solution of the n uncoupled differential equations 


* = (2.4) 


where Q 42 9 are the generalized forces corresponding to qj- 


rjr 

If the functions F,(t) and hence Qj (t) are so irregular that they cannot be 
expressed in mathematical form, equations (2.4) can be solved by Newmark’s 
numerical method. (2) 


The Stiffness and Flexibility Matrices 


Consideration is given here to a multi-story, multi-bay framed structure, 
with girders that are not infinitely rigid, so that joint rotation can take place. 
It is assumed that the structure has a rigid foundation and that the lateral 
movements resulting from the length changes in the columns can be neglected. 

In Fig. 2 a three-story structure of this type is shown schematically. 
Bearing in mind the definition of the stiffness matrix, it is required to find a 
square array of nine numbers, each row of which represents the forces at the 
three levels, compatible with a unit deflection at one level and zero deflection 
at the other two. Moment-distribution is ideally suited to the computation of 
these forces. In Fig. 2 the applied force S99 and the two out-of-balance forces 
Sg1 and S93 form the second row of the stiffness matrix for the structure 


* This and subsequent raised numerals in parentheses refer to the list of 
references at the end of the paper. 
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shown. Only three simple moment distributions have to be carried out for the 
determination of the complete matrix 


$31 $32 $33 


When [S] has been determined and M is known equations of the type (2.2) 
can be written down and solved. This procedure is applicable also to frames 
with semi-rigid connections, provided that for the determination of [S] the 
moment-distribution for semi-rigid frames is used. (3) 

When the effect of the flexibility of the foundation and “beam action” of the 
whole structure—i.e. length changes in the columns—are to be taken into ac- 
count, it is convenient to work in terms of the flexibility matrix rather than 
the stiffness matrix. The former is the reciprocal of the stiffness matrix. 
Stiffness matrices cannot be added by adding algebraically the corresponding 
terms, but this can be done with flexibility matrices. Denoting the flexibility 
matrices due to shear action, beam action of the whole structure and due to 
base rotation by [f'] , [f""], [f''"] , respectively, the total flexibility matrix [f] 
is given by 


Cf] = (ft) + Cen) + [eer] (3.2) 


Some Numerical Examples 


The response of the structure shown in Fig. 3(a) to the forces as shown in 
Fig. 3(b) will now be determined according to several alternative assumptions. 

Although it is indicated in Fig. 3(a) that the moment of inertia of all girders 
is 1675 in.4 the calculations are first carried out on the assumption that the 
girders are infinitely rigid. Moreover, length changes in the columns are 
ignored, as is the flexibility of the foundation. The stiffness matrix for this 
condition is readily computed. The shear K; that must be applied between 
floors to produce a unit relative displacement is 


48EI,, 


K 
r h? 
r 


where I, is the moment of inertia of each of the four columns. Below the 
r-th floor and hy is the distance between floors. Thus 


K, = 166,000 1b. 


K, = K, = 234,000 1b. 


from which 


(s] = (3.1) 

al 

ee 

; 

Wag 
‘ 

Lit 


SCHENKER 
4,00 <-23, 

= }7234 468 -234 14093 


and the equations of motion, after dividing both sides by 1000, are: 
0.832x, + 400x, - 234x, = 52(1 - 0.8t) 


0.832x, + 2 52(1 = 0.&t) (4.2) 
0. 416x5 234x, + 234x3 = 26(1 - 0.8t) 
The solutions of (4.2) are obtainable by the method outlined in section 2. The 


natural circular frequencies are found to be 7.74 sec.~1, 22.56 sec.-1 and 
32.20 sec.~! and the mode shapes in normalized form are: 


1st mode: * 0.0160 0.0241 = 0.0271 
and mode: 912 ~9*9279 = 040027 = 0.0288 


3rd mode: 0.0124 P53 -0.024,8 0.0295 
The uncoupled equations of motion, corresponding to (2.4) are 


59-9a, = 2790(1 - 
G+ 508.5q, = -561(1 - 0.8) 


43 


and the complete solution is 


= 0.0271q, + 0.0288q2 + 0.0295a, 


* 46.5701+1.005 sin(7.7ht - 1.468) - 
= -1.10[1+ sin( 22.56t - 1.536) - 


q3 * 0.12f1+ sin(32.20t -1.546) 0.8t] 


ASCE 944-5 
4 (4.1) 
| (4.3) 
(4.4) 
| + 1036.54, = 122(1 - 0.8t) 
x, = 0,.0160q, - 0.0279q, + 0.012hq3 
x, = 0.0241q, + 0.0027q, - 0.0248q, (4.5) 
“ang 
where 
(4.6) 
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The error introduced by the assumption that the girders are infinitely rigid 
is now examined. It is necessary to compute a revised stiffness matrix al- 
lowing for the fact that the moments of inertia of the girders are as indicated 
in Fig. 3(a). Using moment distribution, the revised stiffness matrix is com- 
puted to be 


372-2270 
fs] «| -227 376 «3751 10° tn.” 


1 (4.7) 


indicating some reduction in stiffness, as is to be expected. From here on, 
the procedure for analyzing the structure is the same as in the previous 
example. Equations of motion, corresponding to (4.2) can be written down 
and solved. The natural frequencies are somewhat lower: 6.57, 19.25, and 
28.73 radians per second. The revised mode shapes are: 


lst mode: = 0.0144 Poy = 0.0241 = 0.0287 
2nd mode: 91,5 = -0.0260 9.. = -0.0038 9.,. = 0.0320 (4.8) 
12 22 32 


3rd mode: 0.0179 -0.0245 933 = 0.0236 


These values should be compared with (4.3). The complete solutions are not 
written down explicitly, but can be obtained as before. 


Behavior in the Inelastic Range 


The dynamic behavior of systems which become inelastic after certain 
critical deflections are exceeded can be treated by methods which are basi- 
cally the same as those used for perfectly elastic systems. A step-by-step 
procedure can be used, re-writing the equations of motion when the stiffness 
matrix changes as a result of yielding in some part of the structure. When 
the stage is reached that unrestrained deformation can take place, a proce- 
dure developed by Salvadori and Bleich(4) may be used. It will be found that 
in this last stage the lowest natural frequency is zero. The appropriate mode 
shapes can be computed as usual. 


The Effect of Damping 


The main purpose of this section is to emphasize the difficulties that are 
encountered in analyzing the response of a multi-degree-of-freedom struc- 
ture, if account is to be taken of the effect of damping. To facilitate mathe- 
matical manipulations it is customary to assume that damping is viscous, i.e. 
that the damping forces are proportional to the rate of change of cisplace- 
ments. For a single-degree-of-freedom system there is then no further dif- 
ficulty. In a multi-degree-of-freedom system a new problem arises: is 
damping a function of the absolute or relative velocities of stations or of 
both? Even if the degree of damping is known or has been assumed—say, 10% 
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of critical damping—what are the values to be used for the damping coeffi- 
cients C in the equations 


= F(t) (6.1) 


The answer, in part, is that in many cases it is not necessary to know the 
vaiues of the coefficients C, provided the percentages of critical damping in 
each mode are known. The procedure then is to solve the equations of mo- 
tion (6.1), neglecting damping, as a first approximation. The effect of damp- 
ing can be introduced subsequently into the uncoupled equations of motion of 
the type (2.4), one of which is associated with each mode. These uncoupled 
equations, with damping, have the form 


As an example, for damping in the first, second and third modes of 10%, 15% 
and 20% respectively, €; = 0.1W , €9 = 0.15W9, €3 = 0.2w3. Having solved 
(6.2) for the qj the actual displacements x, are again found from (2.3). For 
the assumption that the total damping effect is the sum of the damping effects 
associated with the various modes, the mode shapes @¢ are not modified by 
damping. 

In some applications it is desirable, however, to be able to determine the 
numerical values of the coefficients C in equations (6.1). This is the case, 
for example, when the equations are to be solved directly bv means of an ana- 
log computer. The problem, then, is to find the relationship between the coef- 
ficients C—of which there are n“ in number—and the critical damping coeffi- 
cients €;—of which there are n. This may be done as follows. The homo- 
geneous part of the uncoupled equations of motion is 


of which the solution is 


6.4 
a;7e (B, cos wit + Dy sin w jt) (6.4) 


The origin of time can be chosen so that at t = 0 all x, = 0. Then the values 
of X, at this time can be denoted by k,,. From (2.3) and the orthogonality 
condition for the mode shapes 


O, when j # k 
p45? 
1, when j =k 


r=n 
| 
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it is found that 


r=n 
(6.5) 


Therefore the values of qj at t = 0 are zero and all B; in (6.4) are zero also. 
Hence 


xt 
qj = D; e sin w st (6.6) 


The equations of motion with damping can be written 


These must hold at t = 0, when = 9, dijo = Dj; and There- 
fore, at t = 0, (6.7) reads 


{-2beu} + {Dw} = fo} (6.8) 


When (6.8) is expanded it may be seen that each of the n equations contains 
terms associated with all modes. Since it is feasible that an initial displace- 
ment may be imparted to the structure that will excite only one mode, the 
equations resulting from (6.11) must be true for any mode, one at the time. 
Hence 


(MIC = (6.9) 


where [2Dé€w] and [Dw] are diagonal matrices. It is thus permissible to can- 
cel Dw on both sides. Consequently, 


= (cll¢] (6.10) 
where 
| 00 eee O 
0 €,0 O 
° eco (8.11) 
0 00 «ecco En 


(6.10) are n2 equations in which the coefficients C are the only unknowns. In 
expanded form, for n = 3 
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from which C14, Cy9, and Cy can be calculated. There are six more equa- 
tions from which the other C’s can be found. For the numerical example of 
section 4 of this paper, for 10% critical damping in all modes, €; = 0.774 
sec.~1, €, = 2.256 sec.~! and €3 = 3.220 sec.~1, so that 


3433-1210 -198 
[c] 3385 -1151 | sec,~2 
- 198 -1151 1785 


The equations of motion can now be set up in the form of (6.1) and solved with 
the analog computer. 


CONCLUSION 


As outlined in the Introduction, the approach to the solution of analysis and 
design problems involving structures subject to dynamic loads can be broken 
down into several phases. Of these, the present paper deals with methods of 
computing the response of tall structures to lateral dynamic loads. The em- 


phasis is on some aspects of the problem which are often not taken into ac- 2 7 


count, namely (i) joint rotation, (ii) deformation into the inelastic range, and 
(iii) damping. No suggestion is made at this time, as to whether these con- 

ditions will affect practical analysis and design procedures to any large ex- 

tent. The contribution this paper is intended to make is to illustrate how the 
special aspects can be dealt with whenever necessary. 

The procedure of setting up the equations of motion of a multi-degree-of- 
freedom system by making use of the concept of the stiffness matrix is fre- 
quently used in the aircraft structures field in one form or another. It might 
well be adopted by the civil engineer, because it throws a good deal of light on 
the dynamic behavior of a structure whose static properties are known. 

The implications of the usual assumptions concerning damping in a multi- 
degree-of-freedom structure are often not fully appreciated. It is with this in 
mind that the discussion in section 6 of the paper was prepared. 

This paper is based on parts of the author’s doctoral dissertation which 
was developed under guidance of a committee under the chairmanship of 
Bruce G. Johnston, Professor of Structural Engineering, at the University of 
Michigan. Many helpful suggestions were made by Professor M. A. Brull, 
Department of Aeronautical Engineering, University of Michigan. 
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Fig. 2 Force System for the Second Row 
of the Stiffness Matrix of a Three-Story 
Rigid Frame Structure 
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(a) Details of Three- Story Structure 


* 52,000(1 - 0.8t)1bs. 


F, = 26,000(1 - 0,8t)1bs. 


t 
(bo) Varietion of Externgl Forces 


Fig. 3 Deteils of Three-Story Structure and 
Applied Forces 
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THE VISCOUS SUBLAYER ALONG A SMOOTH BOUNDARY 


H. A. Einstein,* M. ASCE, and Huon Li 
(Proc. Paper 945) 


SYNOPSIS 


If the sublayer is visualized as a steady, quasi-laminar flow, great diffi- 
culty is encountered in the physical description of the transition to the turbu- 
lent part of the flow. The proposed model visualizes a periodic growth and 
decay of the sublayer. The magnitude of the life-period of the sublayer can 
be predicted theoretically and has been checked by experiment. As a conse- 
quence of this model some predictions are made on the generation of turbu- 
lence near the boundary. 


INTRODUCTION 


The viscous sublayer in the sense of this derivation is a usually rather 
thin fluid layer along the smooth boundary of a turbulent fluid flow in which 
due to the proximity of the wall the turbulent pattern of flow can not develop 
freely. The presence of this layer under the given conditions is today admit- 
ted by all serious students of turbulence, but its exact physical description is 
not equally well agreed upon. The one point, where general agreement exists 
is the fact that at the boundary the entire shear stress is transmitted, in the 
fluid by viscosity, i.e., by molecular action. Outside the layer most of the 
shear is transmitted by eddy-viscosity or, in other words, by momentum ex- 
change between fluid masses of larger than molecular size. But how the 
transition from molecular to molar units takes effect is not generally agreed 
upon, even if an ever-increasing amount of information becomes available 
from measurements in the proximity of frictional boundaries. 

The engineering description of frictional effects in pipe and channel flow 
becomes most consistent if the assumption is made that turbulence is created 
at the frictional boundary. This assumption is supported by the observation 


Note: Discussion open until September 1, 1956. Paper 945 is part of the copyrighted 
Journal of the Engineering Mechanics Division of the American Society of Civil 
Engineers, Vol. 82, No. EM 2, April, 1956. 

* Associate Prof., Mechanical Eng., Univ. of California, Berkeley, Calif. 

1. Elevated numbers refer to corresponding numbers in the list of literature. 
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that turbulence is created only in connection with a frictional boundary or at 
the interface between fluid masses of different velocity. The latter is com- 
monly found wherever separation occurs of a flow from a solid boundary, 
leaving a fluid body of different velocity between the main flow and the solid 
boundary. A good example for such a case is found along the wake of any sub- 
merged object, where the interface between main flow and wake is not stable, 
but steadily fluctuates, shedding constantly newly created eddies into the flow 
if the Reynolds Number is sufficiently high. The same process may also be 
observed at the wakes of the individual roughness elements of a hydraulically 
rough boundary. With the wakes located between the roughness elements, 
turbulence derived from a rough boundary can thus be said to originate at the 
wall, speaking from the viewpoint of the flow as a whole. 

If the turbulence pattern and, particularly, the velocity distributions are 
compared which are derived from smooth boundaries with those along rough 
boundaries, the patterns are found to be identical, i.e., the equation of the 
time-average velocity distribution near a smooth boundary can be transformed 
into that near a rough boundary by mere replacement of the viscous length 
scale V/u, by a multiple of the roughness k,, (0.300 k,). This shows that 
there exists a basic similarity between the two processes. Particularly, it 
may be expected that also along the smooth boundary the turbulence must be 
created predominantly at the boundary itself, not somewhere in the free flow. 
A satisfactory model of the viscous sublayer must thus also explain how tur- 
bulence is created in this case. It is the purpose of this paper to show that a 
model can be found which describes the flow system along a smooth boundary 
on this basis and which is not contradicted by any available measurements. 


A. The Proposed Model 


To arrive at a workable model, let us imagine first that at a given instant 
the turbulent flow continues all the way to the boundary with a finite main 
velocity at the wall. The high velocity gradient at the wall calls for an ex- 
tremely high viscous shear there which is not transmitted into the free fluid 
due to lack of velocity gradient farther away from the wall, but which decel- 
erates the fluid adjacent to the boundary, creating a viscosity controlled layer 
or, to be consistent with common language, a viscosity controlled sub-layer, 
which grows with time in thickness. As the thickness of the sublayer grows, 
the boundary shear reduces and the growth becomes increasingly slower. 

This sublayer is thus an unsteady laminar or near-laminar flow bounded 
on one side by the rigid wall, on the other side in a rather undetermined man- 
ner by the turbulent flow. It is comparable to a shear flow between two solid 
walls which have a constant velocity difference and which continuously in- 
crease their distance with fluid supplied constantly through the wall represent- 
ing the turbulent fluid at a rate necessary to satisfy the continuity condition. 

It is easily visualized that the Reynolds Number of the sublayer continuously 
grows under these circumstances, indicating an ever-growing tendency to be- 
come unstable. It is to be expected under the steady effect of the disturbances 
from the turbulent flow outside that the sublayer flow will become unstable 

and turn into turbulence itself as soon as a critical Reynolds Number has been 
reached. Then turbulent mixing will occur between the previously turbulent 
flow and the newly created sublayer turbulence. With turbulent mixing much 

‘ more effective than viscosity, it must be expected that the sublayer fluid is 
then accelerated by the turbulent outside fluid to about its own velocity in a 
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time much shorter than that which was required for the viscous build-up of 
the sublayer. After this the cycle is closed and a new cycle begins. 

The inherently unsteady character of this flow is demonstrated by Fig. 1, 
in which two instantaneous views of flow conditions in a plastic transparent 
2" diameter pipe are given. The discharge in both is 0.008 cfs, the average 
flow velocity 0.367 ft/sec., the Reynolds Number 5700, the Weisbach friction 
factor f = 0.036, the friction velocity u, = 0.0246 ft/sec. and the sublayer 
thickness as commonly defined is: 


11.6 
Uy 


A small hole was drilled into the bottom-wall of the pipe at about 6.0 cm of 
the scale permitting the injection of dye into the boundary sublayer. The dye 
was introduced at a very low constant rate and velocity. The dye should be 
expected not to mix with the turbulent flow if the viscous sublayer flow is 
steady. It is seen that in one photograph mixing occurs about 3 cm down- 
stream—in the other almost immediately, indicating the existence of sporadic 
mixing as proposed in the unsteady sublayer model. 


1/2 


Fig. 1. Mixing of dye from the viscous sublayer of turbulent 
pipe flow into the main flow. Upper photograph shows 
almost immediate mixing from tube at cm 6, lower 
photograph shows delayed mixing, depending on phase 
of sublayer. 


The above described cyclic growth and decay of the sublayer represents 
the fundamental idea on which the proposed model is based. Some modifica- 
tions have been introduced, however, for the following mathematical descrip- 
tion. First, the introduction of an imaginary outside boundary which permits 
the passage of fluid is actually not necessary since the distance to which 
viscous effects are felt regulates itself automatically if the depth of the tur- 
bulent flow is assumed to be large compared to the sublayer thickness. The 
displacement thickness of the sublayer may then be introduced as a measure 
for its thickness. 
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Secondly, in order to receive expressions which can be integrated and 
otherwise handled mathematically, the eddy-viscosity of the turbulent flow is 
assumed to be infinitely larger than the molecular viscosity of the fluid during 
the period of decay, resulting in a constant average velocity at all points of 
the turbulent flow. During the period of build-up the effect of eddy-viscosity 
is neglected. 

Thirdly, it is assumed that the area of boundary over which simultaneous 
growth and decay of the sublayer occurs is in all directions large compared 
with the sublayer thickness, reducing the sublayer generation to a one- 
dimensional unsteady flow problem. 

Introducing the following nomenclature, the generation of the sublayer can 
be described on the basis of the given three assumptions: 


x,y,z Cartesian coordinates. 
z=0 smooth boundary, flow at z>0. 
U, V, W average ¥ 2locity in turbulent flow, V = W = 0 at all points. 
u, V, W_ instantaneous sublayer velocities, v = w = 0. 
representative velocity in turbulent range near the sublayer. 
time. 
total period of growth for a sublayer. 
The equation of motion for a thin layer parallel to the wall is then 


du 


(1) 


The boundary conditions for a growing phase of the sublayer beginning at t = 0 
are 


for t = 0 %, for all z>0 


for t>0O and z2=0 u=0 


for t>0O and 


A solution for this flow is 


and h a variable of integration. 
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ely with the parameter H defined as od 
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This solution, which is the Gaussian error integral, actually satisfies the 
boundary conditions (2) as can be seen very easily. For t = 0 the limit H be- 
comes infinite and the integral 


@ 


such that the first condition (2) is fulfilled. The same argument holds for 
t>Oandz=co. For z = 0 andt > 0 the upper limit becomes H = 0 and with 
zero range of the integral of solution (3), the velecity u becomes zero, such 
that all boundary conditions are actually fulfilled. 

The error integral can not be expressed algebraically in closed form, but 
has been tabulated and thus made available for numerical calculations.” 


B. The Thickness Scale of the Sublayer 


Equation (4) already indicates that the instantaneous boundary-sublayer 
thickness should be measured by a sliding scale, i.e., by the unit /Vt, since 
all distances z from the boundary are measured by this unit. Let us try to 
understand the meaning of this unit. For this purpose the so-called “displace- 
ment thickness” of the boundary layer is calculated. The displacement thick- 
ness is the distance by which stream lines far away from the surface are 
diverted due to the velocity reduction of the incompressible fluid between the 
stream line and the boundary. The continuity equation for this flow indicates 
that the displacement thickness §* may be determined as 


dz (6) 
° 


° 


if an equal discharge is assumed to flow at all times between the boundary 
and the stream line. This interpretation is taken from the case of a boundary 
layer which is visualized as long, but not infinitely long, in which case the 
continuity may be established between a cross section 1 ahead of the boundary 
and section 2 through the boundary of Fig. 2. 


Boundary Layer 


Next, the value of 6* is determined, introducing the value of u from eq. 
(3) into eq. (6) 
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@ 


dh) dz (1) 


h:o 


With the help of eq. (4) the distance z is now expressed in H, and eq. (5) is 
used to expand 


Now the order of integrations is changed 


@ 
H=0 


which can once be integrated 


2 
= 4 h dh 
hzo 


A second integration gives 


* vt | 2 


From this may be seen that the distance scale / vt of eq. (4) has physically 
the significance of being /r /2 times the displacement thickness. The re- 
sult that the displacement thickness is proportional to / vt could have been 


derived from dimensional considerations, since /vt is the only parameter 
of the problem with the dimension of a length. But the constant of propor- 
tionality can be obtained only by integration. 


C. The Shear Stress in the Sublayer 


To become more familiar with the character of this flow, the shear stress 
may be calculated. As velocity gradients exist in the z-direction only one 
may write 
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72 
(9) 
rv t 
and at the boundary, where z = 0 
U 
(10) 


At time zero, when the turbulent flow with the constant velocity still reaches 
the boundary, the shear stress is infinite and subsequently reduces fast to 
ever-decreasing values. The significance of this infinite value may be tested 
by calculating its impulse or time integral. It may be agreed that an infinite 
instantaneous T ,, - value may be acceptable in a physical problem as long as 
it has no lasting effect. The time integral is thus determined between time 0 
and a final time T 


T T 


tro te 


if the dynamic viscosity is expressed as ¢ = v .p where p is the fluid density. 
It may very easily be seen from eq. (11) that the impulse for finite time inter- 
vals T is always finite and tends toward zero as T is reduced towards zero, 
indicating the insignificance of the high T, at t = 0. From this the time aver- 
age shear 7, may be calculated as 


=— 
(12) 


and if common nomenclature is used which introduces the friction velocity u, 


as 
(13) 


one can rewrite eq. (12) as 
2 2 
T 


and a relationship between T, u, and U, results in the form 


4 U Vv 
(14) 


4 


+ 
ES 
Po 
a 
Si 
=e 
a 


945-8 EM 2 April, 1956 


This eq. (14) indicates that the growing period of the sublayer may be pre- 
dicted under the assumed conditions from the kinematic viscosity of the fluid, 
the average flow and the shear velocities. 


D. The Decay of the Sublayer 


The decay process has been described already in general terms and the 
time necessary for the decay into turbulence and for the reacceleration of the 
sublayer fluid to the general flow velocity U,, outside the sublayer was as- 
sumed to be short compared with the build-up period T of the sublayer. Any 
values averaged during the build-up period may thus be assumed to represent 
averages for the total time. 

The decay process itself is one of instability and is highly non-linear. At 
this time any mathematical description of the decay process itself appears to 
be utterly impossible. The condition of incipient instability, on the other hand, 
may be studied by a stability condition of the laminar sublayer flow. A dimen- 
sional consideration indicates that this condition is characterized by a critical 
Reynolds Number as in case of the corresponding laminar pipe or channel 
flow. If the Reynolds Number is introduced in the form 


* 


it may be expected to correspond to the Reynolds Number for a pipe flow 
based on the pipe diameter. The use of the displacement thickness instead 
of the hydraulic radius may still change the resulting critical Reynolds Num- 
ber from that of a pipe by a certain amount, and it is quite interesting to do 
some estimating of Np*. This can be done by estimating the velocity Up for 
a given u, from the main turbulent flow. 

One may assume that Up exists at a distance z = né4 from the boundary 
in the main turbulent flow. Then 


Uo 


Ux 


= 5.75 Log,, (9.05 Bree (16) 


If eqs. (14) and (8) are used to express 6 jin eq. (16), Up/u, may be deter- 


mined as 12.39 for n = 1.0 and from that Ne = 780. Similarly, values of NR 
can be obtained for other values of n as shown in Table 1. 


TABLE 1 


CRITICAL REYNOLDS NUMBERS POR VARIOUS ASSUMPTIONS OF U 
| 


U, at y = 0.50 = 525 T=131 ha? 8 
U, aty = 1.08 = 780 T = 195 v/a, 8 15.75 Y/u, 
Uj at y= 2 8 = 1073 T = 268 = 18.49 "/a,, 
Ujaty=3 8 = 1256 T= 31h 20.00 
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It is possible now to express all the variables of the problem in terms of 
this critical Reynolds Number as indicated in Table 1. The choice of the 
proper value for z/ é 7 may be made only empirically from measurements, 


although one may conjecture that Table 1 covers the range within which the 
correct value may be found. 


E. Momentum Transfer Through the Sublayer 


Equations (8) and (12) can be combined into the form 


which may be interpreted as a momentum equation. It indicates that the 
shear stress 7 . which is transmitted from the wall to the sublayer per unit 
wall area is used totally to reduce 1/T times per second the velocity of a 
mass p § T" from U, to zero. An equal amount of momentum is transferred 


according to the same equation by the reacceleration of the same mass to Up 
after it has become turbulent. This explanation and description of the momen- 
tum transfer from the sublayer to the turbulent flow is one of the great ad- 
vantages of this model. It is clear that with this model not only momentum, 
but also fluid, dye, sediment, heat, etc. are transferred from the stream 

lined flow of the sublayer to the main turbulent flow at the end of each period 
T and from the turbulent flow to the sublayer during its build-up. Such trans- 
fer of matter has been observed and can not be explained with a basically 
steady flow pattern in the sublayer. 


F. The Average Velocity Profile in the Sublayer 


If the duration of the decay period is again assumed to be negligible with 
respect to the period of build-up, the average velocity distribution in the sub- 
layer can be determined by averaging the velocity at a given distance z from 
the boundary over the build-up period T only. On the basis of eqs. (3) and (4) 


one can write 
Uo fr’ 


In this form eq. (18) shows that the proposed sublayer has a universally ap- 
plicable velocity distribution if 0/U, is given in terms of z/2 JvT. Herein 
are 


@ the time average sublayer velocity at a distance z from the boundary. 
U, the assumed to be constant turbulent velocity outside the sublayer. 

v the kinematic viscosity of the fluid. 
T the period of the sublayer build-up. 


=. 
(17) 
ag 
| 
it 
ae 
| 
He 
= 


945-10 EM 2 April, 1956 


This function of eq. (18) can not be integrated in closed form. Numerical in- 
tegration gave the curve of Fig. 3 in the form of t/a, against u,z/y. 


G. The Root-Mean-Square (RMS) Deviation of the Instantaneous Sublayer 
Velocity From Its Average 


The RMS value of the local velocity can also be determined numerically 
only. Its determination was based on the following transformation 


T T 
RMS = +/ (u-u) dt = at 
T T 
u* dt-20 at+at fat} 
Tle 


T 
(19) 
T 
° 


This expression can again be made dimensionless by division with u, as is 
done in Fig. 4. 


H. The Rate of Energy Dissipation in the Sublayer 


It is very interesting to determine the percentage of the total flow energy 
which is dissipated into heat in the sublayer by viscous action and the re- 
maining percentage which enters the main flow in the form of newly created 
turbulence. The rate of viscous energy dissipation can be directly calculated 
per unit time and volume 


(20) 


The time average dissipation, again per unit of time and volume, is 


T 
yy du 2 
ree (21) 
if is taken from eq. (9) 
z 
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Uc — TIME-AVERAGE CENTER VELOCITY 


FIG. 4 RMS VALUES OF THE LOCAL FLOW VELOCITY 
IN THE SUBLAYER DETERMINED FOR VARIOUS 
VALUES OF Up AND COMPARED WITH MEAS- 
UREMENTS BY J. LAUFER. 
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2vt (23) 
2 
> ds (24) 


This integral can be numerically solved from tables. The result is given as 
a plot of ‘ 


2 
z 


Particularly interesting is the rate of total energy dissipation in the en- 
tire sublayer. This total rate @ may be obtained by a double integration as 
follows for the unit boundary area 


du 
> vp ( ) dzdt 
~Zz:0 


2t 
t#0 


T 2 
P Uo dt: - Yo yz 25 
PU, us — (25) 


The total rate of energy © t available in the main turbulent flow may be de- 


termined for the unit area and time interms of the average velocity U,, of 
the turbulent flow 


= 2 
= Ugy. = Puy (26) 
from which the fraction of energy dissipation in the sublayer is 


® | Uo (27) 
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I. Application of the Model to Actual Flows 


In its application any such model must describe as many as possible of the 
various phases of the actual flow behavior. This model has the advantage 
over any quasi-steady-flow sublayer description, that is explains the exchange 
of matter between the turbulent flow and the boundary as it has been observed 
to exist (see Fig. 1). This is the great achievement of this approach. 

No model and no theory are better, however, than the assumptions on which 
they are based. While it is believed that the basic description of the sublayer 
cycle is correct, some of the mathematical simplifications used in the deriva- 
tion and solution of the equations are somewhat questionable. The assump- 
tion of an infinitely high eddy-viscosity and of the resulting constant flow 
velocity in the turbulent flow are definitely unrealistic. It will be necessary, 
therefore, to investigate the effect of this idealization of the flow. This can 
best be done by comparing actually observed velocity distributions with the 
prediction of the model. Fig. 3 shows a number of predicted velocity distri- 
butions based on various assumed U, - values and corresponding experiment- 
al values measured by Laufer.? The measurement seems to confirm the 
contention that in a graph (u/u,) against (u,z/y) the distribution is independ- 
ent of the overall Reynolds Number of the flow. The graph indicates, further- 
more, that the outside velocity U, must probably be measured at a distance 
of between 2 and 367 * to give a velocity distribution similar to the measured 
distribution. 

Let us assume for the sake of discussion the value z = 3.06 rT" indicates 


the point of the turbulent velocity distribution curve (shown in Fig. 3 as 
dashed curve) where U, is obtained. At this point u/u, = 15.6, indicated by 
the horizontal tangent to the calculated u/u, - curve. It may be remembered 
now that due to the assumption of uniform turbulent velocity distribution the 
flow velocities in the range of the entire sublayer is assumed to have the 
value u/u, = 15.6 at the beginning of the viscous sublayer build-up period. 
This velocity distribution is assumed to have been created in the short pre- 
ceding period of decay by turbulent momentum exchange with the outside tur- 
bulent flow. It would be more appropriate to assume that at the beginning of 
the viscous sublayer period the velocity distribution follows the extended tur- 
bulent velocity distribution. Unfortunately, this approach leads to hitherto 
unsurmountable mathematical difficulties. Also the physical picture becomes 
rather confused in that case, since then the transmission of momentum is 
simultaneously performed, partially by viscosity, partially by turbulence. All 
we may conclude with security at this time is that the proposed theory as- 
sumes at time zero near the boundary an excessively high flow velocity. With 
the time-average velocity nearly correct, the final minimum velocity is prob- 
ably too low and it must be expected that the predicted RMS-values of the 
velocity near the boundary are too high. 

Fig. 4 shows that this is actually so. Here the curves of the RMS-values 
of the sublayer velocity according to eq. (19) were plotted for various U,- 
values and compared with values measured by Laufer. All curves are high 
along this steep rise from the wall as we have just predicted. They come to 
a sharp maximum near (u, z/y) = 6 - 12 and then drop again reasonably sharp 
towards zero, while the measured values remain high for the rest of the sec- 
tion. This is perfectly reasonable, since the predicted curve describes the 
unsteady laminar flow only, while the measured points include in addition the 
fluctuations of the turbulence. Turbulence must be expected to reach its 
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maximum intensity just outside the sublayer, such that the difference of the 
predicted curves and the measured values is well explained for values of 
(u, z/v) larger than the maximum of the predicted RMS-curves. 

From this comparison of measured and predicted values of velocity mea- 
surements, one must thus conclude that the proposed theory, though basically 
correct, gives only a rough approximation of the actual flow pattern. 

This has a further consequence. The values for the predicted period T 
must also be expected to be larger than the actual values. This may easily 
be seen from the derivation of eq. (17). If the range of velocity change at 
various distances from the wall is effectively less than assumed, the number 
of changes per unit of time must be higher than predicted, to permit an equal 
shear to be transmitted by momentum exchange. Unfortunately, no measure- 
ments could be found in the literature from which the period T could be taken. 
There is not even mention that any periodicity of the process is indicated by 
available measurements. 

The authors of this paper were so convinced of the superiority which the 
unsteady sublayer theory has over any quasi-steady description that they de- 
cided to attempt the experimental proof of its correctness. 


J. The Basis of the Experimental Demonstration of the Periodic Character 
of the Sublayer 


The fact that none of the successful investigators of sublayer flow in the 
past noticed any periodicity is sufficient proof that the periodicity can not be 
directly apparent in any velocity signal, but is probably concealed by large 
amplitude random disturbances originating in the turbulent flow outside the 
sublayer. It appeared thus to be most promising to locate the point of mea- 
surement as far away from the turbulent flow as possible, i.e., to measure 
at the frictional boundary itself. The flow velocity vanishes there and thus 
can not be measured. It must be expected, on the other hand, that the pres- 
sure reflects any such periodicity of the flow pattern. A pressure recording 
system was developed using a Rutishauser pressure gauge as described in 
Appendix I. But also this signal was expected to contain a large percentage 
of random noise created by the turbulence of the main flow, calling for a 
special method to detect any possible periodicity. The use of the available 
sound analyzers was soon dismissed as unreliable since their band width was 
excessive and their frequency range insufficient for low frequencies. Con- 
tinuous records were thus: made of the pressure and an auto-correlation curve 
constructed for the records. This method has the following value: Even if 
the major part of the signal fluctuations are irregular random fluctuations, 
this part of the signal will result in a regularly dropping correlation curve. 
If there is a periodic component added to the random signal, this part of the 
signal will appear as a superimposed periodic fluctuation in the correlation 
curve, too, as may be seen, for instance, in Fig. 6. 

The auto-correlation coefficient for a continuous function of time such as 
the boundary pressure is 


2 
Pit) FP 


where Pt) is the pressure at time t, Pits At) the same pressure at A t later 
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and where the three integrals must cover the same range of t. It is obvious 
that the integral values may be approximated by the sum of a large number 
of individual p-values. It is necessary, however, to choose the number of 
these individual points from 200 to 1000 to eliminate by large numbers the 
accidental influence of individual extreme values. The determination of each 
individual C-value represents from 4 to 12 hours work if this method is used 
without the help of any special tools. Fig. 6 shows that only a large number 
of C-values for a wide variety of A t-values is able to show any characteris- 
tic trends such as the periodic component in Fig. 6. This direct method was 
thus abandoned as too slow. 

According to Putz et al(4,5,6,7) the correlation coefficient C can also be 
obtained in a different way. If a new function 7 is introduced such that 7 has 
the value of +1 wherever p is larger than its time average value, and has the 
value of -1 wherever p is smaller than its time average; then a new correla- 
tion coefficient Cy; may be determined as 


™ (t+ At) at 
C,= 


2 
ot 


Putz (4,5,6,7) shows that Cy = C if p is a random function of the kind of 
most turbulent flow components. This coefficient C; is much more easily 
determined since a simple counting process replaces the much more com- 
plicated multiplication process. But even this method is still quite time- 
consuming —about 2 hours work for an individual Cy-value. But this process 
is very adaptable to mechanization and a device is being developed now which 
is able to produce C,-values in a matter of minutes. Then the analysis of any 
records by correlation curves becomes an easy matter, and much information 
is expected to become available soon on the periodic component of turbulent 
processes near boundaries. The writers had hoped to include some of these 
results in the present paper, but, unfortunately, the development of the me- 
chanical correlation machine turned out to be more time-consuming than an- 
ticipated. Thus only a few curves exist at this time which were derived by 
manual determination, such as those shown in Figs. 6and7. While the periodic 
component is very pronounced in Fig. 6, this tendency is not so clear in Fig. 7. 
The flow of Fig. 6 may thus be used to check the frequency of the fluctuations 
against the theoretical prediction. The flow under consideration is an open 
channel flow of machine oil with the characteristics as given in Appendix I. 
Assuming again that Up must be measured at 3.06 TT the period T may be 


calculated to be T = 1.25 sec. The measured C,-curve shows a distinct 
period of 0.4 sec., or about 1/3 of the predicted value. Whether the above 
discussed idealization of the boundary conditions for the analytic treatment 
can be made responsible for such a large deviation can probably not be de- 
cided until much more experience is available with this method of statistical 
analysis. 


K. Boundary Sublayer and Turbulence 


In order to point out some of the relationships between the sublayer flow 
and the turbulence outside the sublayer, let us begin with the description of 
some of the characteristics of turbulence. Many definitions have been given 
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of turbulence in the past, most of which stress the random character of the 
flow pattern in turbulence. Its detailed description still is made by means of 
the Navier-Stokes equations, however. 

One of the most important characteristics of turbulence can be pointed out 
if a comparison is made between fluid agitated by simultaneous complicated 
wave systems and that agitated by turbulence. The wave systems visualized 
here may be of such amplitude, such that the resulting motion can be de- 
scribed by linearized equations. Again, as in the case of turbulence, a mo- 
tion results which appears to be random in character if a large area is ob- 
served, and which becomes absolutely continuous in the small detail—and 
still, it is not turbulence. The main difference is that in case of the wave 
motion any dye which one may inject into the fluid will for a long time remain 
at about the same location without much mixing with the surrounding fluid, If 
dye is injected in a turbulent flow, however, it mixes very intensely with the 
surrounding fluid and can very soon not be identified any more. This indicates 
that neighboring particles will not stay close together but soon occupy widely 
separated positions. This behavior is usually found in flows with strong shear 
motion or shear deformation of the fluid. Now one can object that a laminar 
flow also shows this deformation, but is neither turbulent nor will it mix dye 
as efficiently with the surrounding fluid as will the turbulent flow. The main 
difference between the laminar and the turbulent flow seems to be that the 
shear motion is regular and even in the laminar flow, but is irregular in the 
turbulent flow. 

This same difference can be expressed also in the following way: Both 
types of motion are of the type in which the fluid has vorticity. Their differ- 
ence is that in the laminar flow the vorticity is distributed regularly accord- 
ing to the main flow pattern, while its distribution in turbulence is in patterns 
which are of much smaller scale than the main flow pattern. High concentra- 
tions of vorticity appear to be characteristic for a vigorous, fresh turbulence. 
There is no mechanism in a turbulent flow which is as well fit to explain the 
origin of these high concentrations of vorticity in the fluid than the proposed 
model of a periodically disintegrating sublayer. 

No explicit expressions shall be given at this time of the vorticity distri- 
bution which periodically originates at the boundary to become part of the 
main turbulence at the end of each period T. It may suffice to say that the 
distribution of vorticity at time T may easily be determined from the equa- 
tions of the sublayer and that this distribution is probably not much changed 
during the process of disintegration of the sublayer. This qualitative discus- 
sion may be sufficient until more data exist to prove or disprove the validity 
of the proposed mathematical description of the unsteady sublayer. 


L. Possible Improvements of the Solution 


It is believed that the given solution is qualitatively correct, but that, 
among others, the following improvements are desirable if possible: 


1) It has already been pointed out that the assumption of a constant velocity 
in the turbulent flow adjacent to the sublayer is very artificial and has been 
introduced only for mathematical reasons. 

2) It was assumed that the entire sublayer becomes part of the turbulent 
flow. It is possible that only the outer part of the sublayer becomes turbulent 
and that an inner part of the sublayer remains always at the boundary. Even 
if this steady part of the sublayer is very thin, it would materially change the 
constants calculated for 6 7*, T, etc. 
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3) It is easily possible that the one-dimensional description of the sublayer 
growth is insufficient to predict the actual conditions if the real growth occurs 
not simultaneously over large areas of the boundary. It is conceivable that the 
process actually changes phase over the various parts of the boundary like a 
wave motion. Different solutions would result in that case. 

4) It is most desirable that a description be found for the disintegration 
process as that part may occupy sufficient time to change materially the time 
average values of average velocity, velocity fluctuation and of the period T. 

An intense qualitative and quantitative study of the processes as well as an 
extension of the proposed theoretical approach will be necessary to find an- 
swers to these questions. 
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SUMMARY 


1) It has been observed that dye and sediment particles are exchanged be- 
tween a smooth boundary and the turbulent flow through the viscous sublayer. 
This fact makes the assumption of a quasi-steady sublayer impossible. 

2) No reasonable explanation can be given for the creation of new turbu- 
lence and of the transmission of shear from a quasi-steady sublayer to the 
turbulent flow. 

3) The periodically growing and disintegrating sublayer pattern is pro- 
posed instead which explains all these facts satisfactorily. 

4) It is assumed that the period of disintegration is negligibly short com- 
pared to the growing period; average values for the entire period are thus de- 
rived from a simplified description of the growing period alone. 

5) The most objectionable assumption made for mathematical reasons is 
that of constant flow velocity in the turbulent flow outside the sublayer. It is 
shown that the deviations of observed patterns from the prediction of the 
theory can qualitatively be explained by this assumption. 

6) No qualitative contradiction has been found to date between the proposed 
theory of the sublayer and flow observations. 
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APPENDIX I 


Experimental Procedure used in the Determination of Local Wall Pressures 


After the problem as such was conceived about three years ago a search 
was made of various pressure recording devices, i.e., devices which can be 
used to record continuously the instantaneous pressure and its fluctuations as 
they apply to a small part of the boundary surface of a turbulent flow. An ap- 
praisal was made of the frequencies which must be expected to occur in the 
problem, of the magnitude of the amplitudes of the fluctuations and of the area 
which might be expected to show predominately equal pressure. On the basis 
of water flowing at Reynolds numbers of about 10* to 105 the following speci- 
fications were drafted for the pressure gage: 


Frequency response 0 to 5000 cps 
Pressure sensitivity 10.001 psi 
Maximum pressure 1 psi 


Area covered by the sensitive element 0.001 Ft. Diameter 


No instrument could be found commercially which satisfied these stringent 
conditions even approximately. It was thus decided to acquire the one which 
approximated the conditions the closest. It was hoped that either the esti- 
mated conditions were too stringent or that the instrument could be improved 
later to adapt it to the particular purpose. 

The best prospect was at the time a Rutishauser transducer. This instru- 
ment uses a membrane of 0.0156 ft. free diameter which has a frequency re- 
sponse up to 10,600 cps and a normal sensitivity of 0.01 psi with a maximum 
pressure of 1 psi using a diaphragm of 0.001 inches thickness. It may easily 
be seen that this instrument is still far from satisfying the original conditions, 
but it appeared to have the best performance characteristics of any instru- 
ments then on the market. The membrane together with an electrode 0.0016 
inches behind it represent a variable condenser which is used to modulate the 
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frequency of a 10 megacycle oscillator. The output of this oscillator is then 
analyzed by a discriminator resulting in a D-C output signal proportional to 
the original motion of the diaphragm. This electric signal was then amplified 
and recorded or its RMS value measured. For the details of the electronic 
system the reader is referred to report (LI, 1954).8 

With this kind of instrument it is possible to adjust the useful ranges some- 
what. If the sensitivity is to be increased, for instance, this can be done by 
replacement of the diaphragm by another one of smaller thickness. Unfor- 
tunately, this change also reduces the frequency response, such that sensi- 
tivity can be gained by increased amplification until the unavoidable tube 
noises become excessive. The area of measurement can be reduced by 
mounting the diaphragm inside a cavity, which itself is connected with the 
free flow by a small opening. Again, the reduction of measuring area is 
traded for frequency response, which is drastically reduced by the cavity. 
For the theory of the cavity, the reader is again referred to report by 1.8 

A very important and extremely time consuming phase of the experimental 
work was the calibration of the pressure indicator and its complex electronic 
amplifying and recording system. This calibration could not be made with one 
piece of equipment for the entire range, but was divided into a number of in- 
dependent stages: 


1) A static calibration was made by direct comparison with manometer 
readings under use of DC-amplifiers. Only limited sensitivity can be 
reached by this method because of the difficulty of stabilizing high gain 
DC amplifiers. 


A low frequency calibration was made by means of a mechanically 
driven excenter which periodically compressed the air enclosed in a 
bellows. Frequencies up to 160 cps could be reached with this instru- 
ment which automatically governed the amplitude independent of the 
frequency. 


A high-frequency vibrator was used for frequencies between 130 and 
6000 cps. This system consists of a loud speaker driving various 
resonance tubes, in the bottom of which the pressure heads were 
mounted. 


The calibration procedure became complicated not only due to the use of 
the different instruments but also because the third method does not permit 
an absolute calibration but only a relative comparison of various instruments. 

The calibration of the transducers with the membranes directly exposed to 
air was performed first. The first two methods which permit an absolute 
calibration showed that both pressure heads with all available diaphragms 
gave a perfectly straight response up to the maximum frequency of 160 cps. 
The theoretical limit of linearity can be predicted to be 3500 cps for the thin- 
nest membrane. It would be expected that if any deviation from the straight 
calibration would occur in the range of frequencies between 160 cps and the 
frequencies where resonance begins to affect the calibration, such deviation 
would change with the instrument and particularly with the thickness of the 
membrane. It could be expected that the calibration was still straight for 
both instruments as long as their signal was the same for a given signal in the 
resonance tube. By means of this method the instruments were shown to have 
a straight frequency response up to a frequency of 2000 cps for all diaphragms 
in air. 
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The next step in the calibration was the introduction of a cavity filled with 
liquid with a reduced opening to the open fluid in which the pressures are 
measured, reducing the area over which the pressure is averaged. The addi- 
tional drastic reduction of the frequency response which the theory of the 
cavity predicts was actually found to occur by this measurement. The table 
on the preceding page gives the frequency response for various geometries 
of the cavity. 


The Experiments 


The previously enumerated values for the qualifications of the measuring 
device are compared with the values as found to characterize the instruments 
by their calibration one finds that 


1) For a given opening of the measuring cavity the frequency response is 
reduced considerably for a significant part of the active spectrum. 

2) The sensitivity of the device is barely sufficient to separate the ex- 
pected signal from unavoidable noise. . 


This state of affairs appeared to prove the chosen instrument to be utterly 
unsatisfactory. Before the entire approach was changed it was decided to try 
if possibly at least certain partial problems or certain ranges of conditions 
could be investigated by the available instrumentation. Thus, the question 
was asked under which conditions large scale of the turbulence pattern near 
the wall is combined with low frequencies and high pressure scales. The only 
possible answer is a flow with low Reynolds number and high viscosity—a flow 
of a highly viscous fluid. After previous attempts had already shown that 
closed pipe systems are very difficult to be freed from all kinds of pressure 
waves which constantly travel through the system and would distort the pres- 
sure records, it was decided to construct an open channel flume circulating 
a flow of oil at high velocity, but moderate Reynolds number. 

The flume with a rectangular cross section of 11.5 inches clear width and 
8 inches depth was constructed of wood, supported by timber bents. The oil 
dropped from the flume into a reservoir from which it was pumped by an 
electrically driven centrifical pump through a 10 inch pipe to the upstream 
end of the flume. The discharge was adjusted by by-passing part of the dis- 
charge directly into the reservoir. One of the most important features of the 
system was the complete mechanical separation of the flume proper from the 
pump and the return pipe to keep all vibrations away from the flume and the 
pressure recording instrument which is very sensitive to vibrations of any 
kind. The discharge was measured at a contraction of the return pipe; the 
oil depth in the flume by point gages. A continuous record was kept of the oil 
temperature, from which the viscosity could be determined after a set of 
viscosimeter tests had determined the viscosity-temperature relationship for 
the particular oil. For engineers who are accustomed to work with water 
flows, the resulting flow was rather odd: with respect to turbulence the flows 
were near critical, i.e., from laminar through critical to turbulent, had about 
2 inches depth of flow and velocities of 6 to 7 ft/sec., and was thus flowing 
with velocities significantly higher than wave velocity, or with other words a 
laminar shooting flow was very readily obtained. 

This experimental arrangement had another great advantage. A laminar 
flow was first established with velocities close to critical and the pressure 
recording instruments turned on such that the noise level could be observed 
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as it was generated by the pump and pipeflow and by other unrelated activities 
in the neighborhood. Then, the flow was slightly increased to cause turbulence 
to be generated and the record continued while the flow measurements were 
made. It could be observed that the pressure records showed fluctuations due 
to turbulence of a larger order of magnitude than the general noise level. 
Experimental Conditions for the Flow of Fig. 6. 
The following flow conditions were established 

Channel width b = 0.957 ft. 

Channel slope 8S, = 0.0502 

Flow depth h =0.144 ft. 

Flow velocity Uy = 6.27 ft/sec 

Discharge = 0.865 cfs 

Viscosity of oil = 1.23 x 1079 ib. sec/ft-2 

Density of fluid =1,73 slugs/ft? 

Kinematic viscosity vy =17.15x 10-4 ft2/sec 


If the shear is assumed to be evenly distributed over bottom and sides, the 
average shear stress T , is 


= 0.313 lb/ft? 
Fo/p = 0.180 ft#/sec” 


Shear velocity u, VT, /p = 0.423 ft/s. 
Thickness of the viscous sub-layer as conventionally defined is 


6 = 11.6 v/u, = 1.96 x 107? ft. 


The period of sub-layer growth and separation frequency are for various 
Uo 

U, at 0.56 ,* T = 0.524 sec n = 1,91 sec 

U,, at 1.0 5 .* T = 0.780 sec n = 1,28 sec 

U, at 2.08,,* sec n = 0.935 sec” 

U, at 3.0 5,.* T=1.25 sec n = 0.800 sec™! 


1 
1 


Since the flow could not be maintained for a very long time for various 
practical reasons it was decided to record continuously the pressure for 
several minutes and to postpone the analysis of the record for a later time. 

The block diagram of the experimental equipment is shown in Fig. 8. 

The Tektronix low level preamplifier type 122, which amplifies the dis- 
criminator voltage, has a noise level as low as 1 to 4 microvolts expressed 
as input signal. The preamplifier has a maximum gain of approximately 1000, 
the frequency response curve is essentially flat from 1 cycle per second to 
40 kilocycles per second. 

The output signal of the preamplifier was further amplified in a brush 
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universal analyzer and recorded by a Brush magnetic oscillograph BL-202. 
The Brush universal analyzer is a D-C amplifier and has a frequency 
response which is flat from zero to 100 cps. For higher frequency signals a 
DuMont Cathode-Ray Oscillograph type 304-A was used to replace the Brush 
universal analyzer, in connection with a DuMont oscillograph recording 
camera type 321. A short excerpt of this record is given in Fig. 5. 
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Block Diagram of Equipment for 
Pressure Fluctuation Measurement of two points. 
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Discussion of 
“VIRTUAL MASS AND ACCELERATION IN FLUIDS” 


by T. E. Stelson and F. T. Mavis 
(Proc. Paper 670) 


T. E. STELSON,! J.M. ASCE, and F. T. MAVIS,2 M. ASC¥.—The writers’ 
data on added mass of rectangular parallepipeds with square side moving 
broadside-on (Fig. 4) have added meaning as they have been replotted by Mr. 
Caldwell. Using the ratio of added mass to displaced-volume mass as 
ordinate he shows quantitatively in his Fig. 5 what every canoeist has 
experienced—namely that for the same acceleration it takes much more force 
to push a paddle through water broadside-on (small relative length) than 
feathered (large relative length). Figs 4 and 5 show the same experimental 
data for rectangular parallepipeds—but Mr. Caldwell’s Fig. 5 shows them 
more clearly. The data can be formulated as follows: 


-0.94 
0.67 é) 


where 


c = added mass 

P = mass density of the fluid 

V = volume of fluid displaced 

d = thickness in the direction of motion 
w = broadside width of square face. 


Mr. Silberman questions these statements in the introduction to the paper: 


“Determination of added mass by experiment often has been inconsistent 
and unreasonable. Measured values for added mass have usually been 
larger than can be explained by theory and analysis.” 


Each sentence is merely a statement of fact as the writers intend it~—and the 
second is not a conclusion drawn from the first as Mr. Silberman may have 
inferred. In fact the data shown in Mr. Silberman’s Fig. 6 eloquently supports 
both statements. 


First, what was measured and plotted in Fig. 6 is expressed as a drag 
coefficient—a variable which mixes velocity-dependent and acceleration- 
dependent forces inseparably. What the writers measured and reported 
in Figs. 2, 3 and 4 of the paper was an “added-force” term that was only 
acceleration-dependent, because in their experiments accelerations were 


1. Asst. Prof. of Civ. Eng., Carnegie Inst. of Technology, Pittsburgh, Pa. 
2. Prof. and Head, Dept. of Civ. Eng., Carnegie Inst. of Technology, Pitts- 
burgh, Pa. 
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high and velocities were low. In fact the writers consider this separation 
and control of variables—variables that have commonly been masked in 
earlier experiments—to be one of their more important contributions to 
experimental studies of virtual mass. 

Second, every point except one as plotted in Mr. Silberman’s Fig. 6 shows 
the “drag coefficient” to be substantially higher than the solid line (1) de- 
signated “Cy, (attributable to added mass). . . .” 


Thus, Mr. Silberman’s Fig. 6 says in the language of graphs the same thing 
as the writers’ two sentences say in the language of words—both essential 
“languages” of the engineer. 

There seems to be abundant other evidence in the literature since 1779 to 
support the statement that “determination of added mass by experiment often 
has been inconsistent and unreasonable.” Yet many early investigators of 
virtual mass deserve great credit for their keen observation and analysis of 
fundamentals that are so elusively hidden by the simultaneous influence of 
uncontrolled variables. One of these was F. W. Bessel* who, in 1828, report- 
ed that because of added mass the period of pendulums was increased from 
one part in 3,000 to one part in 35,000. It is a tribute to Bessel and his ex- 
perimental technique that he was able to detect such slight changes; and yet 
it is not surprising that from such relatively small quantitative measurements 
he should have arrived at deduced added masses that were as much as 90 per 
cent larger than can consistently and reasonably be expected nowadays. 

Incidentally, all but two points shown in Mr. Silberman’s Fig. 6 are from 
100 per cent to 1000 per cent higher than the dotted curve marked “(1) + (2)” 
which, presumably by his hypothesis, should represent the drag coefficient 
for steel spheres accelerated indirectionally. 

The closing suggestion by Mr. Caldwell, namely, to study “. . . the prob- 
lem of a fixed object held in an accelerating flow,”—particularly in a closed 
conduit—may resolve whatever differences of interpretation there may be be- 
tween Mr. Silberman and the writers. It is because of differences in analyti- 
cal and experimental studies that engineers can acquire a better understanding 
of accelerated motion of structures and machines in fluids—and of the associ- 
ated forces they must anticipate in design. 


* Bessel, F. W., “Untersuchungen uber die Lange des einfachen 
Secundenpendels.” Abhandl.d. Akad. d. Wissensch. (Berlin) 1-256 (1828). 
Described in English by F. Baily in Trans. Roy. Soc. London, 399-492 (1832). 
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Discussion of 
“A RESISTOR-NETWORK SOLUTION OF THE ELASTO-PLASTIC 
TORSION PROBLEM” 


by J. H. Weiner, M. G. Salvadori, and V. Paschkis 
(Proc. Paper 671) 


J. H. WEINER,! M. G. SALVADORI,? M. ASCE, and V. PASCHKIS.°—Th> 
writers wish to thank Mr. Moos for his interesting comments and questions 
concerning the paper. His first suggestion for obtaining a more accurate 
value of the roof slope by extrapolation is certainly feasible theoretically. If 
the solutions corresponding to the two mesh sizes were obtained numerically, 
the procedure would be completely practical since it is then known that the 
entire error in each solution (except for the minor round-off error) is due to 
the finite mesh size. However, if the solutions are obtained by a network, 
there is, in addition, a small experimental error which is not a function of 
mesh size, and this error may cause large errors in any extrapolation pro- 
cedure, 

His question regarding the extension of the use of resistor networks into 
other areas is one which has occupied the interest of the Heat and Mass Flow 
Analyzer Laboratory for some time. Development is now under way on their 
use in the solution of thermal stress problems and in the determination of 
stresses in thin shells. The writers hope to be able to report on these ac- 
tivities in the future. 

In regard to his third question, work has continued since the presentation 
of the paper on the direct determination of 2 e; (Eq. 5 ''’), and two success- 
ful techniques have been developed. The measurement of I, - I directly was 
first tried by a circuit, Fig. a, which seems simple and attractive. The re- 
sistance R, equals the equivalent resistance of the network from the voltage 
source, E, to ground if all the nodes of the circuit are grounded. The resist- 
ance Ra is the equivalent resistance of the network with e; at their proper 
values. If rg = rp and if Ra = Ro, the reading of A is 0. Thus by using a 
variable resistor in place of Ro, its proper value can be set by having Ra 
with e; = 0 and adjusting Rp until A = 0. The current in Ra equals that in R, 
and is Ig. Then, without changing the setting of Ro, the values ej and Ra are 
allowed to take their proper values and A is read. The current in Ra is now 
I. One finds - 


4 =m (1) 


1. Asst. Prof. of Civ. Eng., Asst. Technical Director, Heat and Mass Flow 
Analyzer Lab., Columbia Univ., New York, N. Y. 

2. Prof. of Civ. Eng., Columbia Univ., New York, N. Y. 

3. Adjunt Associate Prof. of Mechanical Eng., Technical Director, Heat and 

. Mass Flow Analyzer Lab., Columbia Univ., New York, N. Y. 
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if I, is the current in R, and I the current in Ra. 

However, a mathematical analysis shows that even very small errors, re- 
sulting for example from rg + rp or from A + 0 when adjusting Ro, cause er- 
rors in 4 of the order of magnitude of 20 - 30%. 

Therefore two methods were developed, both of which yield the required 
sum of the nodal values of ej directly. One method consists of an integrating 
device, attached to the voltmeter. In the present case a Leeds and Northrup 
Speedomax recorder, type G, was used through the isolating amplifiers nor- 
mally employed in the Heat and Mass Flow Analyzer. The integrating device 
is mounted on the instrument, which is successively connected to each node 
and permits direct reading of the total voltage, without reading and recording 
the individual values. The integrator actually integrates voltage-time. 
Therefore the meter is kept on each point for a constant length of time, and 
the integral is then divided by this constant time. Although very brief times 
would theoretically be sufficient, it has been found that more reliable results 
are obtained if the constant time on each point is 40 seconds or more. 

The integrator result compared with the sum of voltages obtained by read- 
ing and adding all voltages arithmetically is accurate to within + 1.7%. 

Another method of obtaining voltages directly consists in connecting a very 
high resistance, Ry, to each node. Al! these resistances are then connected 
together and grounded through an additional resistor, Rc, which is calibrated 
accurately (Fig. b). The voltage drop through Rc is proportional to the sum 
of the currents from the several nodes, and each current is proportional to 
the voltage of this node, provided that the voltage at the junction of Rc and 
Ry is negligible in comparison with the nodal voltages ej. Values of Ry 
must be chosen high enough to prevent an influence on the value of e;. It has 


been found that if Ry is 1000 times higher than the resistance R (Fig. 3, 
original paper) results are satisfactory; in the present case Ry, was 10,000 
times higher than Rc. The sum of the voltages thus determined was within 
0.5% of the sum of the individual voltage readings determined arithmetically. 
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Discussion of 
“BLAST RESISTANT BUILDING FRAMES” 


by Bruce G. Johnston and Archie Mathews 
(Proc. Paper 695) 


BRUCE G. JOHNSTON, ! M. ASCE, and ARCHIE MATHEWS,2 ASCE.— 
Mr. McKee has contyibuted interesting suggestions in his discussion and is 
entirely correct in emphasizing that the paper presented an incomplete pic- 
ture of blast resistant frame types in that the specific studies pertain to a 
limited class of buildings. The writers emphasized this limitation themselves 
at the outset but chose this type of building because of its importance to heavy 
industry and because it has not had as much attention as completely protective 
construction. The paper does provide generally applicable suggestions re- 
garding design for maximum structural strength and deformation leading to 
improved plastic behavior of steel beams, steel columns, beam columns, 
compression details, and connections; that are applicable to other types of 
blast resistant steel construction in addition to the completely frangible wall 
type of building that was used as a design example. 

Although, as Mr. McKee points out, the majority of structures have walls 
with windows and may therefore be assumed to be partially open, the writers 
believe that in buildings that are specifically designed for strong blast expo- 
sure it would be desirable either to minimize openings and design the walls 
for complete protection or go as far to the other extreme as is practicable 
and provide walls and/or windows that are frangible. Mr. McKee advocates 
use of shear walls and these are certainly desirable where not prohibited by 
interior accessway requirements. It should also be recognized that shear 
walls will greatly increase the load on the frame if the blast wave strikes the 
structure (with frangible walls) obliquely or fromthe end. From this point of 
view cross bracing would be better and would serve essentially the same pur- 
pose in steel construction as indicated in Figures 4 and 5. Mr. McKee sug- 
gests that the recommended use of reinforced concrete roof slabs may lead to 
collapse under dead load. This is true in some cases for non-continuous 
framing but if fully continuous welded steel construction is adopted, the ex- 
terior columns will always need to be designed for bending as well as direct 
stress since they receive bending moments distributed to them from the ends 
of the roof beams. In such construction it is doubtful if collapse due to dead 
load will ever occur at ten times the yield deflection for the outside columns, 
These deflections would be considerably less than ten times the yield deflec- 
tion for the interior columns, which might be 8WF31 sections as suggested 
by Mr. McKee. 

Mr. McKee suggests the use of members of circular cross section because 
the drag coefficient would be only 0.35. The hollow circular section is an 


1. Prof. of Structural Eng., Univ. of Michigan, Ann Arbor, Mich. x 
2. Research Aide, Dept. of Theoretical and Applied Mechanics, Univ. of Pe 
Illinois, Urbana, Ill. 
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excellent structural member but there are several reasons why it is not used 
to a greater extent. There is difficulty in the design and fabrication of con- 
nections. If welding is used this difficulty may very largely be overcome. 
The greatest deterrent to use of circular members, however, is the fact that 
their cost per pound greatly exceeds that of the open rolled structural shape. 
It is also incorrect to state that the drag coefficient of circular members is 
0.35 without some qualification. The drag coefficient is dependent on the 
Reynolds Number and this is especially true in the case of circular cylindri- 
cal shapes. For example, a 2 1/2" diameter round cylinder with a wind ve- 
locity of 100 miles per hour will have a drag coefficient of 1.20 and it is only 
for a considerable increase in size that the drag coefficient theoretically re- 
duces to 0.35. In practical applications this may be on the low side because 
of the effect of surface roughness in increasing the drag coefficient as com- 
pared to coefficients obtained on machine made wind tunnel models. 

Proposed equivalent static lateral loads for frames supporting frangible 
wall construction were advanced as tentative only. Mr. McKee agrees that 
design for such loads will improve the blast resistance. In an actual war 
emergency the uncertainty as to size, location, and other factors associated 
with possible bomb detonations affecting a building in or very near a target 
area makes the value of precise dynamic analyses questionable. There has 
been no great evidence of desire on the part of builders of industrial plants 
to do very much in the way of strengthening structures against nuclear blast. 
In view of this situation it seems to the writers that best that can be generally 
hoped for at the present time is to introduce those improvements in blast re- 
sistance that can be obtained at a reasonable cost and that will effectively 
reduce the probable area of total destruction. It may be that variations in 
bomb parameters will alter the relative chance for structural survival of a 
specific structure but the over-all result will be to reduce damage in critical 
areas. If the builder can follow the government recommendation to disperse 
his plants to a region at some distance from critical areas it would of course 
be desirable to consider gradually reduced levels cf protective requirements. 
Recommendations for such reduced levels of design load will probably become 
available through future government publications. 

The writers do not belittle in any way the importance of dynamic analyses 
especially as a background for research and for the setting up of standards 
for bomb resistant construction. In the end, however, regardless of what ap- 
proach is used, the engineer will usually select column and beam sizes by use 
of some sort of equivalent static design load. If he uses charts that are 
available as mentioned by Mr. McKee for the solution of simple problems he 
is in effect using the results of prior dynamic analyses to establish the 
equivalent static loads that he will now use to proportion the structure. 

Mr. McKee obviously has a background of study and research in the field 
of blast loaded structures. In general his remarks have been constructive 
and his suggestions as to limitations inherent in some of the ideas expressed 
by the writers are well taken and should receive consideration. 
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Discussion of 
“FLOW INTO A WELL BY ELECTRIC AND MEMBRANE ANALOGY” 


by Chong-Hung Zee, Dean F. Peterson, and Robert O. Bock 
(Proc. Paper 817) 


ROBERT E. GLOVER,! M. ASCE.—The ingenious electric analog arrange- 
ment described by the authors facilitates the investigation of flow conditions 
near a well in an aquifer which has no upper confining bed. Such devices often 
fill a need since the formulas commonly used for engineering purposes are 
based upon simplifications and approximations which exhibit their greatest 
short-comings in this area. 

A common idealization is based upon the following simplifying assump- 
tions: 


a) The gradient of the free surface is conserved throughout the depth of the 
aquifer. 

b) The reduction of the area available for the flow of ground water due to the 
depression of the free surface, as the well is approached, is neglected. 


These two assumptions are equivalent to the single requirement that the draw- 
down everywhere should be small compared to the original depth of water in 
the aquifer. The draw-downs in actual wells often exceed these limitations 
and provide an inducement to seek for more elaborate solutions which will 
provide closer approximations. However, many actual cases involve transi- 
ent phenomena and the added mathematical difficulties introduced by the time 
factor provide powerful inducements for retaining the simple idealization just 
described. The best way out of these difficulties seems to be to use this sim- 
ple idealization as long as the draw-downs retain tolerable values and to sup- 
plement these solutions with a special investigation in the zone of large draw- 
down near the well. The electric analog device described should provide a 
means of ma! ‘ag this special investigation and there is good reason for be- 
lieving that, under certain identifiable conditions, this treatment will be ade- 
quate even though transient conditions must be dealt with. 

If y represents the draw-down, or lowering, of the free water surface in 
an aquifer having an initial water depth D, a permeability K and a ratio of 
drainable voids to total volume v the simple idealization described yields the 
following expression for the condition of continuity: 


For a steady state the right hand member is zero and a comparison with the 
Author’s equation (5) will show that, if the weight of the membrane is neglect- 
ed, the differential equations are of the same form. It follows that the 


1, Research Eng., Bureau of Reclamation, U. S. Dept. of the Interior, Denver, 
Colo. 
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membrane can take forms which represent very closely the draw downs ob- 
tained from this idealization. The steady state solution of equation (a) which 
corresponds to the Author’s equation (1) is: 


where y,, represents the draw down at the well. 

The transient state solution for the case of a constant rate of withdrawal 
Q from a well of radius ry, is given by the well known exponential integral 
formula. 


It has been noted by Taylor and Rainville® that if this integral is expressed in 
the usual series form, in terms of powers of the quantity which appears as 
the lower limit of this integral, plus a logarithm term and a constant, the 
terms of the power series tend to vanish as time increases because the time 
appears in the denominator. The last term to disappear is the one which is 
raised to the first power. If this term is retained and two values of y are 
subtracted an expression of the type of (b) above is obtained together with the 
difference of the two first power terms and it can be seen that the transient 
case reverts to steady state behavior, between the radii r,, and rg, when 


2 
is small compared to 2 (d) 


It is inferred from these results of the approximate theory that if ro 
represents a radius at which the draw-downs remain small compared to D, a 
steady state continuous with the solution (c) at r = rg may be constructed in 
the interval r,, to To, by electric analogy or other means to properly repre- 
sent the true conditions between ry and ro at that time. 

What constitutes satisfaction of the requirement (d) in any given case must 
be decided on the basis of the accuracy desired or attainable. It may be re- 
marked that the aquifers found in nature often depart materially from the 
homogeneity postulated in the idealizations and through these irregularities 
often set a limit to the accuracy which may be attained. 


2. THEIS, C. V. “The Relation between the Lowering of the Piezometric Sur- 
face and the Rate and Duration of Discharge of a Well using Ground Water 
Storage.” Transactions of the American Geophysical Union 1935. 

3. TAYLOR, W. H. and RAINVILLE, E. D. “Analysis of Pumping Test----” 
Informal memorandum., Bureau of Reclamation, Denver, Colorado, 

Feb. 2, 1937. 
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VAUGHN E. HANSEN, ! J.M. ASCE.—It is indeed encouraging to note the 
progress which is being made to solve the complex problem of the hydraulics 
of unconfined wells. Very significant steps have been taken since the writer 
applied the membrane analogy to wells? and introduced the concepts of the 
discharge and well parameters. The experimental data obtained recently by 
several investigators has added considerably to our knowledge of unconfined 
wells. The design curves of figures 9 and 12 are very useful and valuable 
summaries of this experimental data. 

In reviewing the dimensionless plottings of the free-surface profiles shown 
in figure 9, the writer was of the opinion that plots of this nature would be- 
come even more useful and meaningful if parameters could be selected which 
reflected a better picture of the actual free-surface profile. The disadvantage 
of basing the profile upon an arbitrary radius of 115 times the well radius 
could also be avoided by such a procedure. 

In further consideration of this need it appeared reasonable that the follow- 
ing functional relationship of well variables would uniquely define the free- 
surface profile: 


f; (h, Tr, Tw, hg, hy) =0 (1) 


where all symbols are the same as originally defined by the authors. The 
variables of equation (1) can be grouped in dimensionless form as follows: 


r_ tw) =0 (2) 


Since four quantities are still too many to clearly plot on a two-dimension- 
al graph, the following combination of the four parameters of equation (2) was 
considered in an attempt to find three parameters which would give the de- 
sired picture of the free-surface and also provide a unique solution. 


The grouping of parameters shown in equation (3) when plotted with (h-hy)/ry 
as ordinate and r/rw as abscissa with a family of (hg-hw)/rw curves would 
give the desired picture of the free-surface drawdown as the water ap- 
proached the well. The question remained to ascertain if equation (3) would 
yield a unique solution. 

To determine the uniqueness of equation (3), the values of the three dimen- 
sionless parameters were computed for three profiles contained in Mr. Zee’s 


1. Prof. of Irrig. and Drainage Eng., Utah State Agricultural College, Logan, 
Utah. 

2. “Complicated Well Problems Solved by the Membrane Analogy,” by Vaughn 
E. Hansen. Transactions, Am. Geophysical Union, Vol. 33, No. 6, Dec., 
1952. 

3. “Unconfined Ground Water Flow to Multiple Wells,” by Vaughn E. Hansen, 
paper No. 2571, Vol. 118, Transactions, American Society of Civil Engi- 
neers, 1953. 
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Fig. A. Dimensionless Plotting of Free-Surface Profile. 
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Fig. B. Well Parameters of Equation (3) 
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dissertation.4 The results are shown in figure A. It is apparent that two free- 
surface profiles were obtained for the same value of (hg-hy)/rw = 15.0. How- 
ever, when other dissertation data were plotted for r/rw = 115 as shown in 
figure B, a distinct and unique relationship was observed which varifies (with- 
in the range of data available) the validity of equation (3). Unfortunately, the 
three profiles contained in Mr. Zee’s dissertation contained entirely by 
chance two profiles with the same value of (hg-hy)/rw which happened to 
deviate as widely from the best fit line as any of his data. Since additional 
data is not available, the writer is unable to indicate the correct family of 
curves for figure A. However, the proximily and random variation of the data 
from the best-fit curve of figure B does verify that a unique family of curves 
exists for figures A and B. The writer suggests that the authors develop the 
family of curves for figure A with the data used to obtain figure 9. It is the 
writer’s sincere opinion that figure A will be more useful than figure 9 be- 
cause the resulting curves represent the free-surface profile and because 
the arbitrary radius of 115 well radii has been eliminated. 

Attention is called to an error in the fourth line of the last paragraph on 
page 817-10. The symbol h,, should be changed to hg. 


JAMES N. LUTHIN.5—By their ingenious combination of the electric and 
membrane analogies the authors have been able to produce solutions to im- 
portant practical problems of unconfined wells. The method used is simple, 
yields solutions of satisfactory accuracy and is certainly an improvement 
over the laborious numerical methods and other techniques used heretofore. 

A few comments may be in order on two of the boundary conditions which 
are used by the authors. The “free water surface” is an essential part of 
their solution and is provided by the membrane analogy. The “free water 
surface” is normally considered to be in equilibrium with the atmosphere or 
a surface at atmospheric pressure, and is usually assumed to be a limiting 
surface or boundary of the flow region. Muskat(6) defines a free surface as 
a “streamline along which the pressure is uniform.” If we accept Muskat’s 
definition (assuming atmospheric pressure along the free water surface) then 
the occurence of a free water surface in groundwater flow problems is more 
or less fortuitous since streamlines do not necessarily coincide with 
equipressural surfaces. Also the flow system is not limited by the “free 
water surface” since considerable amounts of flow can occur in the capillary 
zone lying above the water table. The water table simply represer’s an 
equipressural surface within the interior of the flow continuum and in general 
is not a stream surface (Luthin and Day, 4). 

The importance of flow in the capillary zone is of course dependent on the 
relationship between the quantity of flow in the capillary zone as related to 
the flow in the groundwater zone or saturated zone lying below the water 
table. Recently Luthin and Day(4) have measured the lateral flow above a 
sloping water table and for the sand used in their experiment where the flow 
section below the water table was of the order of 30 cm thick, about half of 
the total flow was above the water table in the capillary zone. However for a 


4. “The Use of Combined Electrical and Membrane Analogies to Investigate 
Unconfined Flow into Wells,” by Chong Hung Zee, Ph.D. Dissertation, 
Utah State Agricultural College, Logan, Utah, 1952. 

5. Asst. Prof. of Irrig., Univ. of California College of Agriculture, Davis, 
Calif. 
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groundwater zone 50 feet thick the contribution of the capillary zone above it 
would be less than 2 per cent of the total flow indicating that the authors were 
correct in ignoring the capillary zone in their analysis. 

Laplace’s equation is not valid for flow in the capillary region and any an- 
alysis of the total flow continuum depends on a solution of the equation 


div (KV¢) = 0 


where K (the permeability) is a function of the soil moisture tension, The re- 
lationship between K and the soil moisture tension may be determined experi- 
mentally although Childs and Colis-George(1) have proposed a method of cal- 
culating it from the soil moisture characteristic (a plot of soil moisture 
content vs. tension). 

While the assumption of a “free water surface” may be justified in prob- 
lems of well hydrology where the flow region below the water table is very 
large compared to capillary zone lying above the water table, it is doubtful 
if the concept can be applied with accuracy to certain problems of agricultural 
hydrology, highway drainage, and other flow problems involving thin flow re- 
gions below the water table. 

The existence of a seepage surface at the well boundary has been adequate- 
ly treated by the authors as well as others including Hansen(3) but an addi- 
tional comment might be added. As has been pointed out by these authors 
there still seems to be some question in many minds as to the reason for the 
existence of a seepage surface. The usual reason given is that it must exist 
since if we assume the lack of a seepage surface then we have the impossible 
condition of intersecting streamlines with infinite velocities resulting. The 
reasonableness of having a seepage surface becomes apparent when we con- 
sider the condition necessary for water to enter a well. For water to pass as 
a liquid from a porous medium into the atmosphere the only requirement is 
that the pressure of the water in the porous media be equal to or greater than 
atmospheric. [Cf. Richards (71)] Such a condition is most certainly met in 
the situation of ground water flow toward a well. 

In a homogenous soil mass of uniform permeability the length of the seep- 
age surface is independent of the soil permeability. Such uniform soil condi- 
tions are seldom encountered in practice and usually the well penetrates a 
more permeable stratum. The greater the ratio of the permeability of the 
penetrated stratum to the permeability of the overburden the smaller will be 
the seepage surface. This conclusion can be inferred from the results of 
Luthin and Scott(5) who analyzed the steady state flow towards a well pene- 
trating an aquifer having a permeability 100 times that of the overburden. 

The distribution of hydraulic head indicates that the water moves vertically 
out of the overburden to aquifer and then laterally to the well. It seems rea- 
sonable to infer that the seepage surface in this case would be rather small 
since water is passing from a layer of lower permeability to a layer of higher 
permeability [Cf. Day and Luthin (2) ]. 
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CHIA-SHUN YIH.!—The authors are to be congratulated for having con- 
tributed so extensively to well hydraulics. It is, however, desirable to em- 
phasize that the “membrane analogy” really does not exist for unconfined 
wells in the exact sense. For in Eq. (6) of the authors’ paper, the second 
partial derivative of the piezometric head h with respect to the elevation z is 
not equal to zero for unconfined wells. It is true that on the free surface h 
is equal to z, but the equality of h and z is not maintained in the interior of 
the flow, and the error may not be small when the draw-down is not negligible 
compared with the asymptotic depth of the saturated layer. It is precisely 2 
in view of this that the consistency of the authors’ data with those of previous 
investigators—particularly with those obtained by relaxation—is heartening. 
This consistency seems to indicate the general reliability of the authors’ re- 
sults and establishes a posteriori the fact that in spite of the lack of a rigor- 
ous membrane analogy, an aprroximate one nevertheless exists, the use of 
which has made the authors’ extensive data possible. 

Research workers in well hydraulics would render a great service to the 
profession if they would undertake the task of finding the exact solutions for 
a few well chosen cases, with a view to providing additional checks to the 
authors’ results, especially in the regions where only their own data exist. ES 
This would finally and definitely establish the validity of the authors’ exten- ae 
sive results covering a wide range of practical interest. On the other hand, 
the authors could make an additional contribution to well hydraulics if they 
would apply their technique to the treatment of partially penetrating wells or 
of the interference of wells. 


1. Associate Prof., Iowa Inst. of Hydr. Research, Iowa City, Iowa. 
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Discussion of 
“THICK RECTANGULAR PLATES ON AN ELASTIC FOUNDATION” 


by Daniel Frederick 
(Proc. Paper 818) 


HAROLD G. LORSCH, ! A.M. ASCE.—Dr. Frederick’s paper is an excellent 
presentation of an application of the Reissner theory of thick plates. The 
writer fears, however, that the graphs presented might give a distorted pic- 
ture of the numerical errors of the classical theory as compared to the more 
exact Reissner theory used in the paper. 

Figures 4 and 5 show the results for values of (kh4/D) = 1. The author 
states that this value was chosen for convenience. Unfortunately, this value 
does not correspond to usual practical problems. For, using a lean concrete 
plate (E = 2000 ksi, yv = 0.2) and good soil (k = 0.75 kips per cubic inch) the 
thickness of the plate corresponding to (kh4/D) = 1 is almost twenty feet. 
Needless to say that one would expect large shear effects in a plate of that 
thickness. For better concrete and/or softer soil the corresponding plate 
thickness increases. 

For a steel plate (E = 30,000 ksi, v = 0.3) on a thin concrete block 
(E = 2000 ksi, t = 1" = 0", hence k = 200 kips per cubic inch) the plate thick- 
ness for (kh#/D) = 1 must be 13.75 inches, a very unusual value. 

The writer feels that the reduction in bending stress of 55% (line 2, page 
818-12 of the paper) should be viewed in the light of the above numerical cal- 
culations. Plates arising in civil engineering have (kh2/D) ratios from 1/20 
to 1/1000; it would indeed be interesting, if the author presented comparisons 
between the classical and the Reissaer theory for that range. The writer ex- 
pects the two methods to give nearly identical results for those cases. 

None of these remarks are designed to detract from the fundamental ap- 
proach used in the paper. They should merely reassure those engineers who 
might otherwise be inclined to disregard entirely the classical plate theory 
on account of the large errors shown in this paper. 


1. Asst. Prof. of Civ. Eng., The City College, New York, N. Y. 
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